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We investigate the propagation of transverse elastic waves in crumpled media. We set up the 
wave equation for transverse waves on a generic curved, strained surface via a Langrangian formalism 
Oh' and use this to study the scaling behaviour of the dispersion curves near the ridges and on the flat 

facets. This analysis suggests that ridges act as barriers to wave propagation and that modes in a 
certain frequency regime could be trapped in the facets. A simulation study of the wave propagation 
qualitatively supported our analysis and showed interesting effects of the ridges on wave propagation. 
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O ' I. INTRODUCTION 

Waves show distinctive behavior when they propagate through media that are strongly heterogeneous. Among 
these behaviors are localization [Q, multiple propagation paths and lensing . These phenomena depend strongly 
on the type of wave (electromagnetic, Schroedinger, vibrational) and the type of heterogeneity. 

We treat a type of heterogeneity that arises naturally in materials, but whose effect on waves has not been studied: 
crumpled membranes. The disorder introduced by crumpling is much different from the random disorder usually 
considered. Instead of a structureless disorder, crumpled sheets consist of uniform, flat facets bordered by strongly 
curved and stretched borders called elastic ridges |]|l(| . Vibrational waves are expected to propagate much differently 
in the facets and in the border regions. Crumpled sheets of newsprint paper are commonly used as an insulating and 
protective enclosure, but we know of no characterization of their vibrational isolation or transmissions properties. 

We anticipate various possible effects. We imagine oscillating a point on the membrane and investigating the 
response as a function of distance from this source as compared with a flat, uniform sheet. First, there are two 
natural regimes of vibrational frequency to be considered. Low frequencies are expected to produce vibrations with 
wavelengths much longer than the facets or ridges of the crumpled structure. These long-wavelength waves might 
' show distinctive and subtle features, yet we shall not treat them in the present paper. These waves are not expected 
, to depend heavily on the distinctive facet-and-ridge structure of crumpled sheets. 

We shall consider instead the complementary regime of higher frequencies having wavelengths smaller than the 
typical facets of the sheet. Several behaviors might be expected as compared to the uniform sheet. On the one hand, 
the elastic tension induced by the ridges might enhance the propagation along the network of ridges. On the other 
hand, the contrasting nature of the ridges and the facets might serve as a barrier for propagation out of or into a 
facet. It is these phenomena that we address below. 

In section II we briefly review classical thin membrane theory and gather the results we need. In section III we 
derive the wave equations for transverse waves on a generic curved, strained surface from a Langrangian formulation. 
Section IV is devoted to studying the form of the dispersion relations in the two distinct regions (ridges and facets) 
and also the scaling behaviour of the dispersions with thickness of the sheet. This analysis first suggests the possibility 
^ , of trapped energy in the facets with ridges serving as barriers. In section V we discuss in detail the simulations we 
performed and compare them to our analytical predictions. Finally we discuss the implications of our work and the 
scope for future work in section VI. 
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II. THIN MEMBRANE THEORY 



Our findings are based on the classical theory of thin membranes |2|] . If the thickness is small then we may assume 
that all variables may be integrated over the thickness of the membrane and it can be regarded as a two dimensional 
surface ||. To define a surface uniquely up to overall translations and rotations it is sufficient to specify a metric 
tensor g a p and a curvature tensor C a p satisfying Gauss's Theorema Egregium and the Codazzi-Mainardi equations 
H . To see the physical significance of these tensors we first note that the deviation of the metric tensor from the 
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identity is simply the strain tensor 



g a = S a p + -f af 3 (1) 

The sum of the eigenvalues of the strain tensor is the (2D) expansion factor and their difference is the shear angle ||. 
The eigenvalues of the curvature tensor are simply the inverses of the two principal radii of curvature of the surface. 
If we take a point r(xi, X2) where xi, xi refer to the (2D) material co-ordinates of the point, the metric and curvature 
tensor are given by Q 

g a p = (d a r) ■ (dp?) (2) 
C a p = h ■ (d a dpf) (3) 

where d a denotes differentiation with respect to the material co-ordinate x a , and n is the unit normal to the surface 
at that point. We now consider the forces and torques in the system necessary for mechanical equilibrium. They are 
given by the stress tensor a a p and the torque tensor M a p. These are related to the deformation tensors i a p and C a p, 
for sufficiently small strains by || 

Yh r 

C Q /3 = ~ ^ 2 [lap + Vt ap tp T l pT \ (4) 
M a p = K[C a p + VtaptprCpr] (5) 

where Y is the Young's modulus, v is the Poisson ratio and k — Yh? / {12(1 — v 2 )) is the bending modulus of the elastic 
material. Here h is the thickness of the sheet. e a p is the two dimensional antisymmetric tensor. Summation over 
repeated indices is implied. The conditions that guarantee that the strain and curvature tensors describe a physical 
surface in equilibrium are called the von Karman equations. The von Karman equations |^| may be put in a form 
more suitable for our purposes as follows [|| 

det C al3 = dadpiap - V 2 Tr 7 Q/3 (6) 
d a d(3M a( 3 = a a pC a p (7) 

The first equation is simply a restatement of Gauss' theorem which helps in making sure that the curvature and 
metric tensor describe a physical surface. The second von Karman (also called the stress von Karman) equation is a 
condition for physical equilibrium of the sheet. Finally we need expressions for the elastic energy stored in the sheet. 
The stretching energy is given by |(| 

Es = \ J dxdy Gafilafi ( 8 ) 

and the bending energy by 

E b = i J dxdy M a pC a p (9) 
Thus the total elastic energy of the sheet is given by 

E t ot = 2 / dxdy [ a <*plap + M a pC a p] (10) 



III. THE WAVE EQUATION 

Let the equilibrium configuration of the sheet be given by f\)(x,y). This implies that f\)(x,y) satisfies the von 
Karman equations. To study tranverse waves we consider small displacements normal to the surface everywhere. 

f(x,y) = f (x,y) + eu(x,y)h (11) 

where u(x,y) is a smooth function of the material coordinates x,y. Considering only transverse waves is justified 
since it is known that in the thin limit (h <C 1) the transverse and in-plane modes decouple M. Our simulation results 
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also indicate that for small perturbations the generated modes are primarily transverse and the fraction of energy in 
the in-plane modes is orders of magnitude smaller. The transverse perturbation u{x 1 y) leads to a change in the total 
energy of the sheet of order e 2 . The 0(e) term should vanish because fo(x, y) is the equilibrium configuration. Using 
equation (11) along with the defining equations (l)-(5) and equation (10) yields the 0(e 2 ) change in energy 



SE = i J dxdy [k [{u xx - u(C° xx ) 2 ) 2 + [Uyy - u(C° yy ) 2 ) 2 ] + 2/s(l - v)u 2 xy 



hV 

+ 2nu{u xx ~ u(C° xx ) 2 )(u yy - u(C° yy ) 2 ) + y— -^ + u 2 (C° xx ) 2 + v(u 2 y + u 2 (C° yy ) 2 )) 

Hv« + u2 ( C °yy) 2 + « + "'(CD 2 )) + « 2 ((0 2 + Kyf + ^ClA)]] (12) 
where C xx , C y y , 7° x , j y y are simply the respective components of the curvature and strain tensor in the equilibrium 
state. u a is the derivative of u with respect to the material coordinate a. The 0(e) term is simply the stress von 
Karman equation and hence vanishes as argued above. Cross terms involving C xy and ^ xy are absent because we have 
chosen our axes to be aligned with the principal axes of the curvature and strain tensor. Here we restrict ourselves 
to the case where the principal axes of the curvature and strain tensor are the same. This is true along the middle of 
the ridge ]l0|. The general expression is more complex but this is sufficient for our purposes. We now consider SE to 
be an effective potential for the field u{x,y). 

V = SE (13) 



We can then define a kinetic energy term 
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T = - J dxdyph uf (14) 

where p is the density of the material and ut = du/dt. This leads to an effective Langrangian for the field given by 

C = T-V (15) 
The Euler-Lagrange equations for this system are then given by 

dC dC dC 

d^~ d ^ + d ^mr = (16) 

The last term arises because of the presence of second derivatives of the field in the potential Using equations 
(12)-(15) in equation (16) gives us the local wave equation satisfied by the field u(x, y). 

-ph u tt = kV 4 u - A x u xx - A y u yy + Bu (17) 

where 

Yh 

A x = 2k(C° xx ) 2 + 2^(C° yy ) 2 + Y-piJL + n° yy ) (18) 

and similarly for A y . The coefficient of the last term is given by 

B = «(Cl) 4 + K (C° yy r + 2kv{C* x Y{C* v Y + r[f xx {{Cl x ) 2 + „((%,)*) + 7° vy ((C° yy ) 2 + v{C° xx ) 2 )] 

+ r[{Cl x Y + {Cl v f + 2vCl x Cl v ]. (19) 

where r = Yh/ (1 — v 2 ). It is to be noted that all these coefficients are local, in that they depend on the curvature and 
strain at the point (2, y) in the material co-ordinates. We now consider some limiting cases to understand the origin 
of the various terms in equation (17). First we consider the case when we simply have a flat sheet with no curvature 
or strain. Here A x , A y and B are identically zero, giving us 

-ph u tt = kV 4 m (20) 

which correctly describes bending waves on a sheet with thickness h, bending modulus k and density p H. Now we 
consider the case where we have a flat sheet with no curvature, k = 0, and with some isotropic strain 7. Equation 
(17) then reduces to 

-ph u tt = -^!v 2 u (21) 
1 — v 

which describes a stretched membrane under tension proportional to Yj. The leading-order behaviour of the last 
term in equation (17) will be discussed in the next section. 
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IV. DISPERSION RELATIONS 



A crumpled sheet consists of two distinct regions, the ridges and the flat facets that are bounded by ridges. The 
ridges are regions of large curvature and strain and most of the elastic energy of the sheet is stored here H . The flat 
facets on the other hand have small curvature and strain from the ridges pulling on them Jicfl . In this section we look 
at the dispersion relations for tranverse waves in the two regions. 



A. Ridges 

We take the material coordinate x to lie along the width of the ridge and y to lie along its length. We assume that 
the ridge is a region with uniform curvature and strain given by the values that they take at the midpoint. To obtain 
the dispersion relations we consider a plane wave, u ~ e K kx — ut ) j traversing the width of the ridge . Substituting this 
ansatz into equation (17) yields 

afc 4 + bk 2 + c = uj 2 (22) 

(23) 



where 



[ K ^xx + r lxxC 2 x + rv"fy y C 2 x + rC xx + nC yy + 2nvC 2 x C 2 y 



K 

a = ~T 
pn 

b = i- [2kC 2 xx + 2KvC 2 yy + r( lxx + v lyy )\ (24) 
1 

ph 

+ ™ lxx C 2 yy + r lyy C 2 yy + rC 2 y + 2rvC xx C yy \ . (25) 

where r = Yh/(1 — v 2 ) as before. Considering a wave travelling in the perpendicular direction yields the same equation 
except with a x and y interchange in the coefficient b. We now consider the leading order behaviour of the above 
coefficients in the limit of very small thickness , h <C 1, where h is measured in units of the ridge size. For this case the 
scaling of the midridge curvature and strain with thickness have been worked out by Witten and coworkers |9|,|l0|,|6| . 
We have the following results for 

C xx ~ h-i (26) 

Cyy <SC C xx (27) 

Ixx - (28) 

Ivy - ^ (29) 

Using these scaling laws along with the dispersion relation obtained above gives us the scaling behaviour of the 
dispersion curve 

a /i 2 fc 4 + b hik 2 + c h-i = uj 2 (30) 

where ao,6o,co are positive real numbers independent of h. It is also to be noted that we would get the same form 
for the dispersion in the perpendicular direction of propagation. The explicit form of the leading order term in the 
coefficient c is 

Y 



n 2^ ( 31 ) 

where R — C xx is the transverse radius of curvature of the ridge. This is significant because the above expression is 
the square of the "ringing frequency" of a cylinder of radius i?, which is the frequency of purely radial oscillations of 
the cylinder [Q. Thus the last term in the dispersion represents the "breathing mode" of the ridge. Equation (30) 
tells us that, for large wavenumbers (k /i _3 ) ; the first term dominates leading to pure bending wave behaviour. 
This is reasonable, since for small wavelengths (or large k) the surface is locally flat and hence we expect behaviour 
like bending waves on a flat sheet. In the long wavelength regime (fc <ti /i - ^), the constant term dominates leading to 
breathing modes of the ridge. This also prescribes a lower cutoff frequency for modes in the ridge. There is no regime 
where the second term describing sound modes dominates. Thus the crossover leads directly from the bending wave 
regime to the breathing mode regime. 
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FIG. 1. Dispersion curves for the ridge (top curve) and facet (bottom curve) regions for h — 0.01. Coefficients independent 
of h (ao, bo etc) have been set to 1. The lower k cutoff for the ridge region is the inverse of the ridge width (w^ 1 ~ /i -1//3 ). 
The lower cutoff for the facet region is the inverse of the ridge length which we take as unity. Above k w h~ 2 ^ 3 bending modes 
dominate in the ridge and modes in the facet can exist in the ridge too. 



B. Facets 



The facets have some residual strain and curvature due to the ridges. Following the treatment of Lobkovsky et al 
we have the following scaling relations for the curvature and strain away from the ridge. 

(32) 

7~/it (33) 

Using these scaling laws, we follow the same procedure as in the last section to obtain the scaling form of the dispersion 
relation 

ai h 2 k 4 + b 1 hh 2 + Cl h^ =uj 2 (34) 

Again we see two distinct regimes. For small wavelengths (k ^> h~i), we have pure bending waves and for large 
wavelengths (k <C h^z), we have breathing modes. Figure [l] shows the dispersion curves for both the ridges and 
facets. We notice that there are modes that exist on the facets whose frequencies are too low to propagate in the 
ridges. This prompts the question as to whether these modes can tunnel through the ridges to the other facets. 



C. Skin Depth 



To address the question of whether the above mentioned modes are able to penetrate the ridges we look at the skin 
depth of these modes in the ridge region. For a particular mode the frequency is the same in both regions ( because 
of matching at the interface ). If the frequency is less than the cutoff frequency in the ridge region, as is the case with 
the modes we are looking at, then there will be an associated skin depth. This needs to be compared to the width of 
the ridge to see if these modes are trapped in the facet. Consider a mode of frequency luq. In the ridge region we have 

aoh^k 4 + b hik 2 + c h~i = ui 2 (35) 

where ao, Co are just real, positive numbers independent of h. We consider the case where u>q is just less than the 
cutoff frequency. Thus we take 

uj 2 = d h-i (36) 

where d < c o- Using this we may now solve for k in the ridge region. The solution is 

fc= -b ±(b 2 -4a (c -d )) 1 ? ih _, (3?) 
2ao 
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The requirement do < cq ensures that k has a nonzero imaginary part. Thus the imaginary part of k scales as h 3 
and hence the skin depth 

C ~ [Im^)]" 1 ~ h* (38) 

The width of the ridge w, however scales as hs g. Thus we find that 

C < w (39) 

This implies that the modes we considered ought to remain trapped in the facet for sufficiently thin sheets. In the 
next section we examine the predicted scaling behaviour and the trapping phenomenon numerically. 



V. SIMULATIONS 



We model the elastic sheet as a triangular lattice of springs of unstretched length a and spring constant K following 
Seung and Nelson [ fTl) . We incorporate bending rigidity by assigning an energy J(l — rtx ■ ri^) to every pair of adjacent 
lattice triangles with normals ii\ and n^- Sueng and Nelson showed that when the strains are small and the radii of 
curvature are large compared to the lattice constant, this model membrane bends and stretches like an elastic sheet 
of thickness h — ay%J/K with Young's modulus Y — 2Ka/h\/?> and Poisson ratio v — 1/3. The bending modulus is 
k = JV3/2. This model has been extensively used to study the static scaling properties of ridges P, |l0| , p"2"[ . We use it 
in this paper to study the dynamics of wave propagation by adding a unit mass to every lattice point. This models 
the kinetic term in our Lagrangian and assigns to the sheet, a mass per unit area. 

To produce connected ridges and facets we construct a very simple shape, namely the regular tetrahedron, by joining 
the edges of a triangular sheet By using a conjugate gradient minimization routine we get the tetrahedron close 
to its equlibrium position. We then add in damping and a small perturbation (a tap on one of the facets) and then 
wait for the system to attain a state of equilibrium. Once equilibrium is attained we now remove the damping and 
perturb the system ( tapping a few lattice points at the center of a facet by giving the three central particles of the 
facet equal initial velocities normal to the surface) so as to generate waves and study their propagation. To measure 
the energy being carried by the waves we first define an energy that measures the deviation from the equilibrium 
state. 

Edev = ]T X -K{x t - x\f + 2 J(sin(0 fe /2) - sin(0°/2)) 2 

+ E^ 2 (40) 

where x^s refer to the positions of the point unit masses, O^s are supplements of the angles between pairs of plaquets 
and i>i's are velocities of the masses. The superscript indicates the value of the variable at equilibrium. After an 
initial drop (3 — 4%) this energy is conserved up to 1% of its value and hence is a reasonably reliable indicator of the 
energy associated with deviations. 

For large initial perturbations we may also expect the initiation of in-plane waves. To avoid these waves we 
first estimate what initial impulse causes a displacement comparable to the lattice spacing, by assuming that the 
energy of the impulse is transferred to stretching and bending energy in the region of the tap. For typical values 
(a = 1,K = 1, J = 0.05) we obtain the limiting value of the initial impulse to be 0.1. Figure ^ shows the normalized 
energy of the perturbed facet for different tap strengths both above and below the limiting value. The simulation 
was done with the above mentioned parameter values for a tetrahedron of size 10a. We notice that the curves below 
the threshold lie on top of each other. Thus tapping below the threshold strength gives us robust, characteristic 
curves that are independent of the initial impulse. From here on we consider only simulations where we tap below the 
threshold. Looking at the stretching and bending energies in this regime also tells us that stretching contribution is 
small and hence the modes are mostly transverse. We also check this fact by looking explicitly at the particle motion. 

Now we consider the time dependence of the energy of the perturbed facet. Typical curves are like those in figure |^ 
with initial impulse below the threshold. For an initial interval the energy remains close to 1. It then begins to drop 
off and we see a straight line segment to the curve. The initial interval simply indicates the time it takes the fastest 
mode to reach the edge of the facet. The straight line segment on the other hand is presumably where all modes are 
contributing giving us a constant energy flux out of the facet and hence a constant slope in the energy versus time 
graph. 
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FIG. 2. Normalized energy of the perturbed facet as a function of time. E p stands for the deviation energy of the perturbed 
facet. E t ot is the total deviation energy of the tetrahedron. Curves are for values of impulse 0.01 (+), 0.001 (*) and 0.3 (o). 
Note that the first two are below the threshold (0.1) and the last one is above. 
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FIG. 3. Log-log plot of onset times as a function of thickness for different p. Circles (p = 0.9), Pluses (p — 0.85), Asterisks 
(p = 0.92). The lines are best fit to the data points. Slopes of the lines (top to bottom) are -1.55,-1.48 and -1.41 respectively. 

To measure these two quantities we first define an onset time to where the energy begins leaking out to the other 
facets. We take this as the time for the energy of the perturbed facet to drop to a fraction p ~ 0.9. We choose 
0.9 because it is close enough to 1 so as to justify the name onset time but still far enough from 1 to get a good 
spread of data. We now wish to study how this onset time scales with the thickness of the sheet. We obtain data 
for different thicknesses by running the simulation with the parameter a kept fixed and changing J, K such that the 
Young's modulus remains constant while the bending modulus obeys k ~ h? . Note that since the unit masses are not 
changed the density per unit area (ph) remains constant. The onset time should be proportional to the inverse of the 
largest group velocity. From figure |l| we see that the fastest modes should be predominantly bending modes. Hence 
we expect 

h ~ i^r 1 - (v^M)- 1 - h-^ (4i) 

Figure |^ shows a plot of the onset times for various thicknesses and for different values of p. Best fit lines for slightly 
different values of p are fairly parallel and hence the exact choice of p does not affect the power law crucially. From 
the graph we deduce a power law of to ~ jj— i.48±o.07 w } 1 j c i 1 agrees very well with the theory. The error quoted is the 
variation of the slope with p. 

The scaling derived above does not reflect the presence of ridges. The onset time is a characteristic of the velocity 
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FIG. 4. Log-log plot of slope m as a function of thickness. The line is a best fit to the data points. 




10 20 30 40 50 60 



Time 

FIG. 5. Normalized energy of the perturbed facet as a function of time. E p stands for the deviation energy of the perturbed 
facet. Etot is the total deviation energy. Solid line - tetrahedron. Dashed line - flat sheet. 

of propagation in the facet and we would get an identical scaling law for a flat sheet. To identify the effects of the 
ridge we consider the portion of the curve that shows the energy leaving the facet. As mentioned before there is a 
straight line segment to this curve where we may assume all modes are contributing. We denote the slope of this line 
segment as m and study how this quantity varies with thickness. We also expect that if the slopes (m) are dominated 
by the bending wave contribution then 

m ~[^]~^/ph~h?' 2 (42) 

Figure |J shows the expected behaviour with a power h 1A6±01 . Again this is exactly what we expect for a flat sheet. 
This could imply that the ridges do not affect any of the characteristics of wave propagation. 

To see if this really is the case we explicitly compare curves obtained for a tetrahedron and a flat sheet under 
appropriate strain in figure |{| The curve for the flat sheet is obtained by using a flat triangular sheet obtained by 
unfolding the tetrahedron. The sheet is thus divided into four equal "facets" with the central one being perturbed. 
The sheet is put under strain by stretching the edges slightly and then clamping them. The magnitude of strain is set 
to be equal to that in the facets of the tetrahedron. The curves show the normalized energy of the same perturbed 
facet in the tetrahedral and flat geometry. From figure ^| we can draw one main conclusion. The rate at which energy 
leaves the facet seems appreciably less in the case of the tetrahedron. 

We can imagine three possible scenarios which inhibit the flow of energy from the perturbed facet. Firstly, significant 
reflections may occur at the ridge thus confining energy to the perturbed facet. Secondly, the energy may preferentially 
propagate through the centers of the ridges where the curvature is least, i.e. the large curvatures close to the vertices 
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FIG. 6. Normalized energy of a section of the perturbed facet as a function of time. E p stands for the deviation energy of 
the subsection of the perturbed facet. E to t is the total deviation energy. Solid line - tetrahedron. Dashed line - flat sheet. 

may inhibit propagation. Thirdly the energy may simply propagate preferentially along the ridge and not into the 
neighbouring facet. The first two mechanisms follow qualitatively from our analytical treatment while the third would 
invalidate our analysis. To check for the possibility of reflections we look at the energy of a smaller triangular section 
within our perturbed facet. Figure || shows the comparison between the tetrahedral case and the flat sheet case. 
We immediately notice a bump in the curve for the tetrahedron indicating that energy flows back in, suggesting the 
occurence of reflections. 

However to note the presence or absence of other features we looked explicitly at the spatial energy distribution as 
a function of time. Figure (^) shows a sequence of snapshots of the spatial energy distribution for the tetrahedron 
on the left and a flat sheet on the right. The tetrahedron has been folded out and hence the central triangular face, 
marked at the vertices by arrows, is the perturbed face. Ridges form the sides of this face. Arrows in the flat sheet 
sequence simply help to identify the same region. Note that the initial excitation is confined to the central three 
particles. We see that in the flat sheet case, the energy simply spreads outwards. In the tetrahedral case though, we 
notice reflections (t = 6.5 to t = 9) and also the preferential propagation of energy through the middle of the ridges 
(t = 15.5 ). Thus we see the qualitative effects predicted by our analysis. However we never see the complete trapping 
of energy. This could be due to the limited size of the simulation as well as the difficulty in producing wavelengths in 
the right regime. 

VI. CONCLUSION 

In this paper we have taken the first step toward understanding the dynamics of crumpled elastic sheets. We find 
that the ridge structures in crumpled sheets qualitatively modify the propagation of transverse waves. Despite the 
stress induced by the ridges, our analysis indicates that the stretched string or drumhead modes normally associated 
with a stretched membrane are not important in a crumpled sheet. Our analysis suggested that the ridges would 
act as barriers to propagation and also that there could be modes that are trapped in the facets. These modes 
are unable to penetrate the ridges because the associated skin depth is much less than the ridge width in the thin 
limit. Our numerical simulations gave us qualitative evidence that the above picture was right. Future work could 
focus on subjecting our predictions to a more thorough simulation analysis which could support our findings in a 
more quantitative fashion as well as identify the trapping of modes in the right frequency regime. It would also be 
interesting to extend the analysis to cover the low frequency regime. 
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FIG. 7. Time sequence of energy density profiles for the tetrahedron (left) and flat sheet (right). For purposes of display 
the lattice points of the tetrahedron have been mapped to their undeformed triangular state. Arrows mark the vertices of the 
perturbed facet of the tetrahedron and the same region on the flat sheet for comparison. The brightest regions in each snapshot 
represent the regions of highest energy density in that snapshot. 
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